
library("rstudioapi")     
setwd(dirname(getActiveDocumentContext()$path))
getwd()
library("readxl") 
require(miceadds)
library(data.table)
library(tidyverse)
library(DescTools)
library(htmlTable)
library(stargazer)
library(xlsx)
library(estimatr)

load("FinalFinal20062023Full1-6.RData")
load("FinalFinalSP20062023Full1-6.RData")

library(dplyr)
final_finalV2 = mutate(final_finalV2, 
                       Agec = (Age - 30)/10, # per 10 year age
                       Incomec = (Income - 100)/50, # per 50 income
                       clinic_distancec = (clinic_distance - 5)/5, # per 5
                       SNMetricc = SNMetric/10)


# Models 9 with the new data: 
L_9c <- glm.cluster(ActVacApril ~ individual_treatment  + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123)
L_9c2_6 <- glm.cluster(ActVacApril ~ individual_treatment  + Gender + Agec + 
                       Education + Incomec + Employed + SNMetricc + clinic_distancec+
                       dist2 + dist3 + dist4 + dist5 , 
                     data = final_finalV2[which(final_finalV2$`District number`!= "6"),], family=binomial,
                     cluster = final_finalV2[which(final_finalV2$`District number`!= "6"),]$Q123)
L_9c2_5 <- glm.cluster(ActVacApril ~ individual_treatment  + Gender + Agec + 
                         Education + Incomec + Employed + SNMetricc + clinic_distancec+
                         dist2 + dist3 + dist4 +  dist6, 
                       data = final_finalV2[which(final_finalV2$`District number`!= "5"),], family=binomial,
                       cluster = final_finalV2[which(final_finalV2$`District number`!= "5"),]$Q123)
L_9c2_4 <- glm.cluster(ActVacApril ~ individual_treatment  + Gender + Agec + 
                         Education + Incomec + Employed + SNMetricc + clinic_distancec+
                         dist2 + dist3 +  dist5 + dist6, 
                       data = final_finalV2[which(final_finalV2$`District number`!= "4"),], family=binomial,
                       cluster = final_finalV2[which(final_finalV2$`District number`!= "4"),]$Q123)
L_9c2_3 <- glm.cluster(ActVacApril ~ individual_treatment  + Gender + Agec + 
                         Education + Incomec + Employed + SNMetricc + clinic_distancec+
                         dist2 +  dist4 + dist5 + dist6, 
                       data = final_finalV2[which(final_finalV2$`District number`!= "3"),], family=binomial,
                       cluster = final_finalV2[which(final_finalV2$`District number`!= "3"),]$Q123)
L_9c2_2 <- glm.cluster(ActVacApril ~ individual_treatment  + Gender + Agec + 
                         Education + Incomec + Employed + SNMetricc + clinic_distancec+
                          dist3 + dist4 + dist5 + dist6, 
                       data = final_finalV2[which(final_finalV2$`District number`!= "2"),], family=binomial,
                       cluster = final_finalV2[which(final_finalV2$`District number`!= "2"),]$Q123)
L_9c2_1 <- glm.cluster(ActVacApril ~ individual_treatment  + Gender + Agec + 
                         Education + Incomec + Employed + SNMetricc + clinic_distancec+
                          dist3 + dist4 + dist5 + dist6, 
                       data = final_finalV2[which(final_finalV2$`District number`!= "1"),], family=binomial,
                       cluster = final_finalV2[which(final_finalV2$`District number`!= "1"),]$Q123)

tab_9c2_1 = exp(cbind("Odds ratio" = coef(L_9c2_1), confint.default(L_9c2_1, level = 0.95)))
tab_9c2_2 = exp(cbind("Odds ratio" = coef(L_9c2_2), confint.default(L_9c2_2, level = 0.95)))
tab_9c2_3 = exp(cbind("Odds ratio" = coef(L_9c2_3), confint.default(L_9c2_3, level = 0.95)))
tab_9c2_4 = exp(cbind("Odds ratio" = coef(L_9c2_4), confint.default(L_9c2_4, level = 0.95)))
tab_9c2_5 = exp(cbind("Odds ratio" = coef(L_9c2_5), confint.default(L_9c2_5, level = 0.95)))
tab_9c2_6 = exp(cbind("Odds ratio" = coef(L_9c2_6), confint.default(L_9c2_6, level = 0.95)))

rownames(tab_9c2_1) = c("Intercept", "Health", "High Cash", "Low Cash", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access", "District 3", "District 4", "District 5","District 6")
rownames(tab_9c2_2) = c("Intercept", "Health", "High Cash", "Low Cash", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access", "District 3", "District 4", "District 5","District 6")
rownames(tab_9c2_3) = c("Intercept", "Health", "High Cash", "Low Cash", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access","District 2", "District 4", "District 5","District 6")
rownames(tab_9c2_4) = c("Intercept", "Health", "High Cash", "Low Cash", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access","District 2", "District 3", "District 5","District 6")
rownames(tab_9c2_5) = c("Intercept", "Health", "High Cash", "Low Cash", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access","District 2", "District 3", "District 4","District 6")
rownames(tab_9c2_6) = c("Intercept", "Health", "High Cash", "Low Cash", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access","District 2", "District 3", "District 4", "District 5")

# convert function: 
convert_fun = function(data, name_mod){
  new_vec = c()
  for (i in c(1:dim(data)[1])){
    new_vec[i] = as.character(paste(format(round(data[i,1],2), nsmall = 2),paste("(",paste(format(round(data[i,2],2), nsmall = 2), format(round(data[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
  }
  new_vec = data.frame(new_vec)
  rownames(new_vec) = rownames(data)
  colnames(new_vec) = name_mod
  return(new_vec)
}

col_9c2_1 = convert_fun(tab_9c2_1, "Model1")
col_9c2_2 = convert_fun(tab_9c2_2, "Model2")
col_9c2_3 = convert_fun(tab_9c2_3, "Model3")
col_9c2_4 = convert_fun(tab_9c2_4, "Model4")
col_9c2_5 = convert_fun(tab_9c2_5, "Model5")
col_9c2_6 = convert_fun(tab_9c2_6, "Model6")

tb1 = merge(col_9c2_1, col_9c2_2, by = 0, all = TRUE);rownames(tb1) = tb1[,1];tb1[,1] = NULL
tb2 = merge(col_9c2_3, col_9c2_4, by = 0, all = TRUE);rownames(tb2) = tb2[,1];tb2[,1] = NULL
tb3 = merge(col_9c2_5, col_9c2_6, by = 0, all = TRUE);rownames(tb3) = tb3[,1];tb3[,1] = NULL
ttb1 = merge(tb1, tb2, by = 0, all = TRUE);rownames(ttb1) = ttb1[,1];ttb1[,1] = NULL
ttb2 = merge(ttb1, tb3, by = 0, all = TRUE);rownames(ttb2) = ttb2[,1];ttb2[,1] = NULL

tb_exp = ttb2[c("High Cash", "Low Cash","Health", "Access" ,"Age", "Male", "High-Educate", "Medium-Educate", "Low-Educate","Employed","Mean-Food","Social Media","District 2", "District 3", "District 4", "District 5","District 6", "Intercept"),]
stargazer(tb_exp , type = "text", summary = FALSE)
stargazer(tb_exp , type = "latex", summary = FALSE)

no.obs = c((L_9c2_1$glm_res$df.null), (L_9c2_2$glm_res$df.null), (L_9c2_3$glm_res$df.null), (L_9c2_4$glm_res$df.null), (L_9c2_5$glm_res$df.null), (L_9c2_6$glm_res$df.null))
aic_mod = round(c(L_9c2_1$glm_res$aic, L_9c2_2$glm_res$aic, L_9c2_3$glm_res$aic, L_9c2_4$glm_res$aic, L_9c2_5$glm_res$aic, L_9c2_6$glm_res$aic),2)
log_lik = round(c(as.vector(logLik(L_9c2_1$glm_res)), as.vector(logLik(L_9c2_2$glm_res)), as.vector(logLik(L_9c2_3$glm_res)), as.vector(logLik(L_9c2_4$glm_res)), as.vector(logLik(L_9c2_5$glm_res)), as.vector(logLik(L_9c2_6$glm_res))),2)
mod_spec = data.frame(rbind(no.obs, aic_mod, log_lik))
rownames(mod_spec) = c("Observations", "Akaike Inf. Crit.", "Log Likelihood")
colnames(mod_spec) = colnames(tb_exp)


# ExtData Table s4: 
tb_exp2 = rbind(tb_exp, mod_spec)
stargazer(tb_exp2 , type = "text", summary = FALSE, out="SM_TableS4.txt")
stargazer(tb_exp2 , type = "latex", summary = FALSE, out="SM_TableS4.txt")


